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The development of analytic-gradient methodology for excited states within conventional time- 
dependent density-functional theory (TDDFT) would seem to offer a relatively inexpensive alterna- 
tive to better established quantum-chemical approaches for the modeling of photochemical reactions. 
However, even though TDDFT is formally exact, practical calculations involve the use of approxi- 
mate functionals, in particular the TDDFT adiabatic approximation, whose use in photochemical 
applications must be further validated. Here, we investigate the prototypical case of the symmetric 
CC ring opening of oxirane. We demonstrate by direct comparison with the results of high-quality 
quantum Monte Carlo calculations that, far from being an approximation on TDDFT, the Tamm- 
Dancoff approximation (TDA) is a practical necessity for avoiding triplet instabilities and singlet 
near instabilities, thus helping maintain energetically reasonable excited-state potential energy sur- 
faces during bond breaking. Other difhculties one would encounter in modeling oxirane photody- 
namics are pointed out but none of these is likely to prevent a qualitatively correct TDDFT/TDA 
description of photochemistry in this prototypical molecule. 



I. INTRODUCTION 

A complete understanding of photochemistry often 
requires photodynamics calculations since photochemi- 
cal reactions may have excess energy and do not fol- 
low the lowest energy pathway, and since dynamic, as 
much as energetic, considerations often govern how a re- 
action jumps from one electronic potential energy sur- 
face (PES) to another. Unfortunately, full photodynam- 
ics calculations are prohibitively expensive for all but the 
smallest molecules unless simplifying approximations are 
made. One such approximation is to adopt a TuUy-type 
mixed quantum /classical surface-hopping (SH) trajec- 
tory approaclP21 where the electrons are described quan- 
tum mechanically and the nuclei classically. Then, since 
the quantum part of the calculation remains the compu- 
tational bottleneck, further simplifications are required 
for an efficient electronic structure computation of the 
photochemical dynamics. 

Time-dependent density-functional theory (TDDFT)'^ 
represents a promising and relatively inexpensive 
alternative to accurate but more costly quantum chem- 
ical approaches for on-the-fly calculations of excitation 
energies. TDDFT methods for excited-state dynam- 
ics are not yet a standard part of the computational 
chemistry repertory but are coming on line. In partic- 
ular, the basic methodology for the computation of an- 
alytic gradi ents has n ow been implemented in a num- 
ber of codej^*^ ^ l ^ l ^ l ^° l and important practical progress 



is being made on t he problem of calculating SH inte- 
grals within TDDFltIIl^' i3 | i4 | i5 | i6 | indeed first applica- 
tions of mixed TDDFT/classical SH dynamics are being 
reported"'^'^ Problems in this new technology can be 
expected and modifications in the basic TDDFT methods 
will be necessary to enlarge the class of its potential ap- 
plications. In particular, practical TDDFT calculations 
employ approximate functionals whose performance in 
describing photochemical problems requires careful vali- 
dation. 

Few applications of TDDFT to problems rele- 
van t to photoche mistry can be found in the litera- 
l-^^ ji7 | i8 | i9 | 20 | 2i | 22 | ^^le few attempts made to as- 
sess the value of TDDFT for photoch emical appli- 
cations re port either encouragin ^^'^ l ^^ l ^^l or discourag- 
ing resulta^SEZMMMl xhe fundamental reasons for 
these descrepencies stem from the fact that conven- 
tional TDDFT has several problems which may or may 
not be fatal for modeling a given photochemical re- 
action. The main difficulties encountered in conven- 
tional TDDFT include the underestimation of the ion- 
ization threshold^'', the underestimation of charge trans- 
fer excitations'^^ and the lack of explicit two- and 
higher-electron excitation^^HMMHTl^ Still other difficul- 
ties are discussed in pertinant review j^ ^^P^^^^I'^^T'^^P'^l 

To help TDDFT become a reliable tool for modeling 
photochemistry, we must first obtain a clear idea of what 
are its most severe problems. It is the objective of this 
article to identify the most critical points where improve- 
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ment needs to be made, at least for one prototypical 
molecule and one type of reaction path. The reaction 
chosen for our study is the photochemical ring opening 
of oxirane. The photochemistry of oxirane and of oxirane 
derivatives has been much studied both theoretically and 
experimentally (Appendix |A|. It is a text-book molecule 
used to discuss orbital control of stereochemistry in the 
context of electrocyclic ring opening. Textbook discus- 
sions usually focus on CC ring opening (see e.g. Ref.f^pp. 
258-260) but a discussion of CO ring opening may also 
be found {Rei?^pp. 408-411). At first glance, this would 
seem to be a particularly good test case for TDDFT. 
That is, oxirane is a simple molecule with a relatively 
simple photochemistry where TDDFT should work and, 
if not, where its failures will be particularly easy to ana- 
lyze. In this paper, we will focus on the symmetric ring 
opening of this molecule even though the photochemi- 
cal ring opening does not follow a symmetric pathway. 
The reason for our choice is the usual reasorPS, namely 
that the use of symmetry greatly facihtates analysis and 
hence the construction of and comparison with highly 
accurate quantum Monte Carlo results. A mixed quan- 
tum/classical SH trajectory stud y o f asymmetric ring 
opening will be reported elsewhere^. 

In the next section, we give a brief review of DFT and 
of TDDFT. In sec. |III| we review the formalism behind 
the more exact theory against which we will be comparing 
our TDDFT results. Computational details are given in 
Sec ^ 



IV 



Sec. 



VI 



Section |V| reports our results and discussion and 
summarizes. 



where is the HF exchange energy and E^'^^ and 

E^^^ are generalized gradient approximation (GGA) ex- 
change (x) and exchange-correlation (xc) energies. The 
coefficient Cx controls the amount of Hartree-Fock ex- 
change, being unity for Hartrce-Fock^ zero for pure DFT, 
and fractional (typically around 0.2333) for hybrid func- 
tional. For more info rmation about DFT, we refer the 
reader to Refs.l^SESHni. 

Time-dependent density-functional theory offers a rig- 
orous approach to calculating excitation energies. Un- 
like time-dependent Hartree-Fock (TDHF), TDDFT is 
formally exact. Consequently, although approximate 
exchange-correlation functionals must be used in prac- 
tice, we can hope that good approximations will lead 
to better results than those obtained from TDHF which 
completely lacks correlation effects. In practice, the 
use of hybrid functionals means that TDDFT contains 
TDHF as a particular choice of functional (i.e., the 
exchange-only hybrid functional with = 1). 

The most common implementation of TDDFT is via a 
Kohn-Sham formalism using the so-called adiabatic ap- 
proximation which assumes that the self-consistent field 
responds instantaneously and without memory to any 
temporal change in the charge density. This allows the 
time-dependent exchange-correlation action quantity in 
TDDFT to be replaced with the more familiar exchange- 
correlation energy from conventional TDDFT. Excita- 
tions may be obtained from the linear response formu- 
lation of TDDFT (LR-TDDFT). A key quantity in LR- 
TDDFT is the exchange-correlation kernel, 



II. (TIME-DEPENDENT) 
DENSITY-FUNCTIONAL THEORY 

The (TD)DFT PES for the Ith excited state is cal- 
culated at each nuclear configuration by adding the Ith 
TDDFT excitation energy to the DFT ground state en- 
ergy, so Ei = Eq + LUi. [Hartree atomic units (Ti = m = 
e = 1) are used throughout the paper.] Even though 
the use of parentheses in the expression (TD)DFT bet- 
ter emphasizes the hybrid DFT -I- TDDFT nature of the 
calculation, we will usually follow the common practice 
of simply refering to TDDFT PESs. The purpose of this 
section is to review those aspects of TDDFT most nec- 
essary for understanding the rest of the paper. 

The hybrid DFT + TDDFT nature of (TD)DFT cal- 
culations indicates that problems with the ground state 
PES can easily become problems for the excited-state 
PESs. Thus it is important to begin with a few words 
about the ground-state problem. The simplest methods 
for treating the ground state are the Hartree-Fock (HF) 
method and the Kohn-Sham formulation of DFT. In re- 
cent years, the use of hybrid functionals has permitted 
the two energy expressions to be written in the same 
well-known form. 
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which, along with the Hartree kernel, f'^^{r,r') = 
1/ |r — r'l, and the Hartree-Fock exchange kernel (whose 
integrals can be written in terms of the Hartree ker- 
nel) determines the linear response of the Kohn-Sham 
self-consistent field in_the adiabatic approximation. In 
Casida's formulation'^'*, the excitation energies are ob- 
tained by solving a random phase approximation (RPA)- 
like pseudo-eigenvalue equation. 



(2.3) 



where the matrices A and B are defined by, 
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and the integrals are in MuUiken charge cloud notation, 
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It is to be emphasized that the TDDFT adia- 
batic app roxim ation includes only dressed one-electron 
excitation^2lIlZI This is particularly easy to see in the 
context of the Tamm-DancofF approximation (TDA) to 
Eq. (2.3 1 whose TDDFT variant is usually attributed to 
Hirata and Head-GordoiP^. The TDA simply consists 
of neglecting the B matrices to obtain, AXj — cuiXj. 
The number of possible solutions to this equation is the 
dimensionality of A which is just the number of single 
excitations. In fact, the LR-TDHF TDA (exchange-only 
hybrid functional with Ca; = 1) is simply the well-known 
configuration interaction singles (CIS) method. 

An issue of great importance in the context of the 
present work is triplet instabilities. These have been first 
analyzed by Bauernschmitt and Ahlrichs^^ in the context 
of DFT, but the explicit associati on wi th LR-TDDFT ex- 
citation energies was made lateif^^ESl Following Ref.l^, 
we suppose that the ground-state DFT calculation has 
been performed using a same-orbitals-for-different-spin 
(SODS) ansatz and we now wish to test to see if releas- 
ing the SODS restriction to give a different-orbit als-for- 
different-spin (DODS) solution will lower the energy. To 
do so we consider an arbitary unitary tranformation of 
the orbitals, 



i/.^^(r) = exp [iA + ^/)] V^^(r) 



(2.6) 



where R and / are real operators. The corresponding 
energy expression is. 



i?^(A-B)i?+/f(A + B)/ 



0(A3), 
(2.7) 

where matrix elements of the R and / operators have 
been arranged in column vectors and the 0{X) term dis- 
appears because the energy has already been minimized 
before considering symmetry-breaking. The presence of 
the terms (A ± B) shows the connection with the pseu- 
doeigenvalue problem (2.3 1 which can be rewritten as the 
eigenvalue problem. 



(A + B) (A - B) Zi = Lu]Zi 



(2.8) 



where Zj — Xi — Yj. It is an easy consequence of 



Eq. (2.4) that the matrix (A — B) is always positive def- 
inite for pure functionals (cx = 0), as long as the aufbau 
principle is obeyed. However (A 4- B) may have negative 
eigenvalues. In that case, the energy E\ will fall below Eq 
for some value of /, indicating a lower symmetry-broken 
solution. At the same time, this will correspond to a 
negative value of toj (i.e., an imaginary value of lu].) A 
further analysis shows the corresponding eigenvector cor- 
responds to a triplet excitation, hence the term "triplet 
instability" for describing this phenomenon. In principle, 
when hybrid functionals are used, singlet instabilities are 
also possible. 

Thus imaginary excitation energies in LR-TDDFT are 
an indication that something is wrong in the description 



of the ground state whose time-dependent response is be- 
ing used to obtain those excitation energies. As pointed 
out in Refs.l^ ancP^, the TDA actually acts to decouple 
the excited-state problem from the ground-state problem 
so that TDA excitations may actually be better than full 
LR-TDDFT ones. Typically, full and TDA LR-TDDFT 
excitation energies are close near a molecule's equilibrium 
geometry, where a single-determinantal wave function is 
a reasonable first approximation and begin to differ as 
the single-determinantal description breaks down. 

Interestingly, there is a very simple argument against 
symmetry breaking in exact DFT for molecules with a 
nondegenerate ground state. Since the exact ground- 
state wave function of these molecules is a singlet be- 
longing to the totally symmetric representation and so 
with the same symmetry as the molecule, the spin- up and 
spin-down charge densities will be equal and also have 
the same symmetry as the molecule. It follows that the 
same must hold for the spin-up and spin-down compo- 
nents of the exact exchange-correlation potential. Since 
the potentials are the same, then the molecular orbitals 
must also be the same for different spins. Thus, no sym- 
metry breaking is expected in exact DFT for molecules 
with a nondegenerate ground state. Note however that 
the functional is still dependent on spin so that the spin 
kernels fl^{r,r') and /J^^(r,r') are different. 

The interested reader will find additional information 
about TDDFT in the recent book of Ref.El. 



III. QUANTUM MONTE CARLO 

The quahty of our (TD)DFT calculations is judged 
against highly accurate quantum Monte Carlo (QMC) 
results. For each state of interest, the QMC results are 
obtained in a three step procedure. First, a conventional 
complete active space (CAS) self-consistent field (SCF) 
calculation is performed. The resultant CASSCF wave 
function is then reoptimized in the presence of a Jas- 
trow factor to include dynamical correlation and used in 
a variational Monte Carlo (VMC) calculation. Finally, 
the VMC result is further improved via diffusion Monte 
Carlo (DMC). 

To properly describe a reaction involving bond break- 
ing, it is necessary to adopt a multi-dcterminantal de- 
scription of the wave function. In the CASSCF method, 
a set of active orbitals is selected, whose occupancy is 
allowed to vary, while all other orbitals are fixed as ei- 
ther doubly occupied or unoccupied. In a CASSCF(n,m) 
calculation, n electrons are distributed among an ac- 
tive space of m orbitals and all possible resulting space- 
and spin-symmetry-adapted configuration state func- 
tions (CSFs) are constructed. The final CASSCF(n,m) 
wave function consists of a linear combination of these 
CSFs, like in a full configuration interaction (CI) calcu- 
lation for n electrons in m orbitals, except that also the 
orbitals are now optimized to minimize the total energy. 
In the present article, we consider 4 electrons in an active 
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space of 6 orbitals which represents the minimal active 
space for a proper description of all 9 states of interest. 

When several states of the same symmetry are re- 
quested, there is a danger in optimizing the higher states 
that their energy is lowered enough to approach and mix 
with lower states, thus giving an unbalanced description 
of excitation energies. A well-established solution to this 
problem is the use of a state averaged (SA) CASSCF ap- 
proach where the weighted average of t he ene rgies of the 
states under consideration is optimizecP^'^. The wave 
functions of the different states depend on their individ- 
ual sets of CI coefficients using a common set of orbitals. 
Orthogonality is ensured via the CI coefficients and a 
generalized variational theorem applies. Obviously, the 
SA-CASSCF energy of the lowest state will be higher 
than the CASSCF energy obtained without SA. In the 
present work, we need to apply the SA procedure only 
for the 2^Ai state in the C2v ring opening. Although 
the CASSCF and SA-CASSCF l^Ai energies are found 
to be really very close, we calculate the 2^Ai energy as 
EcAsil'Ai) + [E^isi^^Ai) - A similar 

procedure is used in the case of the corresponding QMC 
calculations. 

QMC method^^ offer an efficient alternative to con- 
ventional highly-correlated ab initio methods to go be- 
yond CASSCF and provide an accurate description of 
both dynamical and static electronic correlation. The 
key ingredient which determines the quality of a QMC 
calculation is the many-body trial wave function which, 
in the present work, is chosen of the Jastrow-Slater type 
with the particular form, 

^Y''''^^?^'Y[j{n„nA,r,A), (3.1) 

where rij denotes the distance between electrons i and j, 
and TiA the distance of electron i from nucleus A. We use 
here a Jastrow factor J' which correlates pairs of electrons 
and each electron separately with a nucleus, and employ 
different Jastrow factors to describe the correlation with 
different atom types. The determinantal component con- 
sists of a CAS expansion, ^fp^^, and includes all possible 
CSFs obtained by placing 4 electrons in the active space 
of 6 orbitals as explained above. 

All parameters in both the Jastrow and the determi- 
nantal component of the wave function are optimized by 
minimizing the energy. Since the optimal orbitals and ex- 
pansion coefficients in ^fp^^ may differ from the CASSCF 
values obtained in the absence of the Jastrow factor J', 
it is important to reoptimize them in the presence of the 
Jastrow component. 

For a wave function corresponding to the lowest state 
of a given symmetry, we follow the energy-minimization 
approach of Ref.^^. If the excited state is not the lowest in 
its symmetry, we obtain the Jastrow and orbitals param- 
eters which minimize the average energy over the state of 
interest and the lower states, while the linear coefficients 
in the CSF expansion ensure that orthogonality is pre- 
served among the states^®. Therefore, the wave functions 



resulting from the state-average optimization will share 
the same Jastrow parameters and the same set of orbitals 
but have different linear coefficients. This scheme repre- 
sents a generalization of the approach of Ref.'^ where 
only the orbitals were optimized and orthogonality was 
only approximately preserved. The present approach is 
therefore superior to the one of Ref.^^, which was how- 
ever already giving excellent results when tested on sev- 
eral singlet states of ethylene and a series of prototypical 
photosensitive moleculea^^'^. We note that, when a CAS 
expansion is used in the absence of the Jastrow compo- 
nent, the method is analogous to the CASSCF technique 
for the lowest state of a given symmetry, and to a SA- 
CASSCF approach if the excited state is not the lowest 
in its symmetry. 

The trial wave function is then used in diffusion Monte 
Carlo (DMC) , which produces the best energy within the 
fixed- node approximation [i.e., the lowest-energy state 
with the same zeros (nodes) as the trial wave function]. 
All QMC results presented are from DMC calculations. 



IV. COMPUTATIONAL DETAILS 
A. SCF and TD Details 

Calculations were carried out with two different com- 
puter programs, namely Gaussian 03^ and a develop- 
ment version ofDEMoN2KP. We use the two programs 
in a complementary fashion as they differ in ways which 
are important for this work. In particular, Gaussian 
can carry out HF calculations while deMon2k has no 
Hartree-Fock exchange. Similarly, Gaussian can carry 
out time-dependent HF (TDHF) and configuration in- 
teraction singles (CIS) calculations which are not imple- 
mented in deMon2k. In contrast, the Tamm-Dancoff 
approximation (TDA) for TDDFT may only be per- 
formed with deMon2k. The present version of de- 
Mon2k is limited to the local density approximation for 
the exchange-correlation TDDFT kernel, while the kernel 
is more general in Gaussian. Finally, GAUSSIAN makes 
automatic use of symmetry and prints out the irreducible 
representation for each molecular orbital. The particu- 
lar version of DeMon2k used in the present work does 
not have this feature so we rely on comparison with the 
output of Gaussian calculations for this aspect of our 
analysis. 

All Gaussian results reported here are calc ulated 
with the extensive 6-311-H-HG**(2d,2p) basis seiP^ll. 
All calculations used default convergence criteria. The 
DFT calculations used the default grid for the exchange- 
correlation integrals. 

The same basis set is used for deMon2k calculations 
as in the Gaussian calculations. The GEN-A3* density- 
fitting basis set is used and density-fitting is carried 
out without imposing the charge conservation constraint. 
The SCF convergence cutoff is set at 10~^. We always 
use the FIXED FINE option for the grid. Our implemen- 
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tation of TDDFT in DeMon2k is described in Ref.l^. 

We use two different density functionals, the local den- 
sity approximation (LDA) in the parameterization of 
Vosko, Wilk, and Nusair-- (referred to as SVWN5 in 
the Gaussian input), and the popular B3LYP hybrid 
functionaP^. The corresponding TDDFT calculations 
are referred to as TDLDA and TDB3LYP respectively. 



B. CAS Details 

All CASSCF calculations are performed with the pro- 
gram GAMESS(USj^. In all SA-CASSCF calculations, 
equal weights are employed for the two states. 

We use scalar-rclativistic energy-consistent Hartree- 
Fock pseudopotentials^^ where the carbon and oxygen 
Is electrons are replaced by a non-singular s-non-local 
pseudopotential and the hydrogen potential is softened 
by removing the Coulomb divergence. We employ the 
Gaussian basis set^^ll constructed for these pseudopo- 
tentials and augment them with diffuse functions. All 
calculations are performed with the cc-TZV contracted 
{llsllpld) /[3s3pld] basis for carbon and oxygen, aug- 
mented with two additional diffuse s and p functions with 
exponents 0.04402 and 0.03569 for carbon, and 0.07376 
and 0.05974 for oxygen. The d polarization functions 
for carbon and oxygen are taken from the cc-DZV set. 
For hydrogen, the cc-DZV contracted (10s9p)/[2slp] ba- 
sis is augmented with one s diffuse function with expo- 
nent 0.02974. 



C. QMC Details 

The program package CHAMF^is used for the QMC 
calculations. We employ the same pseudopotentials and 
basis sets as in the CASSCF calculations (see Sec. IVB I. 

Different Jastrow factors are used to describe the cor- 
relation with a hydrogen, an oxygen and a carbon atom. 
For each atom type, the Jastrow factor consists of an 
exponential of the sum of two fifth-order polynomials of 
the electron-nuclear and the electron-electron distances, 
respectiveljil^. The parameters in the Jastrow factor and 
in the determinantal component of the wave function are 
simultaneously optimized by energy minimization follow- 
ing the scheme of Ref.'^, where we employ the simple 
choice ^ = 1. For the excited states with the same sym- 
metry as the ground state, the ground- and excited-state 
wave functions are optimized in a state-average manner 
with equal weights for both stated^. An imaginary time 
step of 0.075 H^^ is used in the DMC calculations. 
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V. RESULTS 

The Woodward-Hoffmann (WH) rules for orbital con- 
trol of symmetry in electrocyclic reactions were an im- 
portant motivation in the 1970s and 1980s to seek 
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FIG. 1; a) Woodward- Hoffmann orbital correlation scheme. 
Symmetry labels are for the C2v point group in the case of 
reactants and products, for the Cs point group along the dis- 
rotatory pathway, and for the C'2 point group along the conro- 
tatory pathway, b) Thermal ring opening, c) Photochemical 
ring opening. Symmetry labels depend upon the labels of the 
X, y, and z coordinates. In this article, the COC ring lies in 
the {y, z)-plane and the z-axis coincides with the C2 axis, in 
agreement with the lUPAC conventiorP^. 



stereos pecific photoch emical ring-opening reactions of 
oxirane^ ^^ l ^'^ l ^^ l ^^ l ^^ l ^^ l Since some oxiranes (notably 
diphenyl oxirane) do appear to follow the WH rules for 
thermal and photochemical ring opening, the WH rules 
might seem like the obvious place to begin our study 
of symmetric ring opening pathways in oxirane. This is 
somewhat counterbalanced by the fact that oxirane it- 
self is an exception to the WH rules for photochemical 
ring opening and by the fact that it is now clear that the 
WH rules do not apply nearly as well to photochemical 
as to thermal reactions (Appendix [A| . Nevertheless, we 
will take the WH model as a first approximation for un- 
derstanding and begin by describing the model for the 
particular case of oxirane. 
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Many models at the tim e that W oodward and Hoff- 
mann developed their theorj-' ^-'^^--'^ were based upon sim- 
ple Hiickel-like 7r-electron models, so, not surprisingly, 
the historical WH model for oxirane uses a three-orbital 
model consisting of one p-orbital on each of the oxy- 
gen and two carbon atoms. (See Fig. [l]) Elementary 
chemical reasoning predicts that the closed cycle should 
form three molecular orbitals, a, n, cr*, in increasing 
order of energy. Similarly the open structure has the 
"particle-in-a-box" orbitals familiar from simple Hiickel 
theory. Woodward and Hoffmann observed that a reflec- 
tion plane (cr) of symmetry is preserved along the dis- 
rotatory reaction pathway while a C2 rotation symme- 
try element is preserved along the conrotatory reaction 
pathway. This observation allows the reactant-product 
molecular orbital correlation diagram to be completed by 
using the fact that the symmetry representation of each 
orbital is preserved within the relevant symmetry group 
of the molecule along each reaction path (WH princi- 
ple of conservation of orbital symmetry) and connecting 
lowest orbitals with lowest orbitals. Figure [l] shows that 
a conrotatory (con) thermal reaction connects ground- 
state configurations in the reactant and product while a 
disrotatory {dis) thermal reaction connects the reactant 
ground-state configuration with an electronically excited- 
state configuration. Thus, the con thermal reaction is ex- 
pected to be prefered over the dis reaction. Also shown in 
Fig. [1] is the photochemical reaction beginning with the 
n ^ a* excited state. In this case, the dis mechanism 
is expected to be prefered since the con mechanism leads 
to a still higher level of excitation in the product than 
in the reactant molecule while the level of excitation is 
preserved along the dis pathway. 

Let us now see whether this orbital model is refiected 
first in the vertical absorption spectrum of oxirane, then 
in the C2V ring opening pathway, and finally in the con 
and dis ring-opening pathways. At each step, the per- 
formance of TDDFT will be assessed against experiment 
and against other levels of theory. 



A. Absorption Spectrum 

The first step towards calculating the electronic ab- 
sorption spectrum of oxirane is the optimization of the 
geometry of the gas phase molecule. This was carried 
out using the HF method and DFT using the LDA and 
B3LYP functionals. The calculated results are compared 
in Fig. [2] with the known experimental values. It is seen 
that electron correlation included in DFT shortens bond 
lengths, bringing the DFT optimized geometries into con- 
siderably better agreement with the experimental geome- 
tries than are the HF optimized geometries. This bet- 
ter agreement between DFT and experiment also holds 
for bond angles. There is not much difference between 
the LDA and B3LYP geometries with the exception of 
the COC bond angle which is somewhat better described 
with the B3LYP than with the LDA functional. 
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FIG. 2: Comparison of HF/6-311G"(2d,2p), B3LYP/6- 
311G**(2d,2p), and LDA/6-311G'* (2d,2p) optimized geome- 
tries with the experimental gas phase geometry from ReLl^. 
Note that the structure has C2v symmetry. 



TABLE I: Principal oxirane singlet excitation energies (eV) 
and oscillator strengths (unitless). 

TDHF TDLDA TDB3LYP EbqjL 

9.14 (0.0007) 6.01 (0.0309) 6.69 (0.0266) 7.24(s) " " " 
9.26 (0.0050) 6.73 (0.0048) 7.14 (0.0060) 7.45(w)'' 
9.36 (0.0635) 6.78 (0.0252) 7.36 (0.0218) 7.88(s)", 7.89(s) ' 
9.56 (0.0635) 7.61 (0.0035) 7.85 (0.0052) 

9.90 (0.0478) 7.78 (0.0304) 8.37 (0.0505) 
9.93 (0.0935) 8.13 (0.0014) 8.39 (0.0168) 

8.15 (0.0405) 8.40 (0.0419) 

"Gas phase UV absorption spectrurrPSl 
'Obtained by a photoelectric techniqud^^. 
''Gas phase UV absorption spectrurJ^. 



Several experimental studie s of the abs orption spec- 
trum of oxirane are 

availablePlSilSSISfilfil], 

The posi- 
tions of the principal electronic excitations were identified 
early on but their assignment took a bit longer. The ac- 
cepted interpretation is based on the study by Basch et 
Q]pS\ who combined information from vacuum ultraviolet 
spectra, photoelectron spectra, and quantum chemical 
computations, and assigned the observed transitions to 
0(n) 3s and 3p Rydberg transitions. Let us see how 
this is refiected in our calculations. 

While theoretical absorption spectra are often shifted 
with respect to experimental absorption spectra, we ex- 
pect to find two strong absorptions in the low energy 
spectrum, separated from each other by about 0.65 eV. 
Vertical absorption spectra are calculated using TDHF, 
TDLDA, and TDB3LYP, all at the B3LYP-optimized ge- 
ometry. The results are shown in Table |T] TDB3LYP is 
the method in best agreement with experiments, yield- 
ing two strong absorptions red-shifted from experiment 
by about 0.5 eV and separated by 0.67 eV. The TDLDA 
method is qualitatively similar, apart from a stronger 
red shift. In particular, the TDLDA method shows two 
strong absorptions red-shifted from experiment by about 
1.1 eV and separated by 0.78 eV. In contrast, the TDHF 
spectrum is blue-shifted by about 2 eV and is otherwise of 
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180° 




FIG. 3: B3LYP MOs. Left: ring structure. Right: open 
structure. 
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FIG. 4: Walsh diagram for C21, ring opening calculated at the 
B3LYP level. To construct this diagram, the COG bond angle 
was varied and all other geometric parameters were relaxed 
within the constraint of C2v symmetry. The HOMO is the 
26i orbital on the left hand side and the 4^2 orbital on the 
right hand side. 



questionable value for interpreting the experimental spec- 
trum. In what follows, we have decided to assign the low- 
est three absorptions using the results of our TDB3LYP 
calculations. 

This first involves an examination of the B3LYP molec- 
ular orbitals (MOs). These orbitals are shown in Fig. |3] 
The electronic configuration of the ring structure is, 

• • • [6ai(a)]2[26i(n)]2[7ai(3s)]"[462(a*)]° • • • . (5.1) 

The group theoretic MO labels (oi, 02, bi, and 62) cor- 
respond to representations of the C'2v symmetry group. 
The additional labels, cr, n, and a*, show our chemical 
interpretation of the B3LYP orbitals and their correspon- 
dence with the MOs in the WH three-orbital model. 

Confirmation of the chemical nature of our B3LYP 
MOs was obtained by constructing a Walsh diagram for 
C2V ring opening. This graph of MO energies as a func- 



FIG. 5: B3LYP MOs imphcated in the principal UV absorp- 
tions. 



tion of ring-opening angle (a) is shown in Fig. 4] As 
the ring opens, the C{2p) 6ai((T) bond breaks and so in- 
creases in energy. At the same time, the C{2p) 462(0-*) 
antibond becomes less antibonding and so decreases in 
energy. The 0{2p) lone pair 2bi{n) is not involved in 
bonding and so maintains a roughly constant energy 
throughout the ring-opening process. Experimental in- 
formation about occupied MOs is available from electron 
momentum spectroscopy (EMS) via the target Kohn- 
Sham approximatioiP^. The ordering of the B3LYP oc- 
cupied orbitals is consistent with the results of a recent 
EMS study of oxiran^^. In particular, the two high- 
est energy occupied orbitals are seen to have a dominant 
p-type character. The interpretation of the unoccupied 
orbitals is more problematic in that the B3LYP unoc- 
cupied orbitals in our calculations are not bound (i.e., 
have positive orbital energies). We are thus attempt- 
ing to describe a continuum with a finite basis set. It 
is thus far from obvious that the unoccupied orbital en- 
ergies will converge to anything meaningful as the finite 
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basis set becomes increasingly complete. The Walsh dia- 
gram shows that the 4^2 (c*) unoccupied orbital becomes 
bound for ring opening angles beyond about 80°. It is 
also rather localized and hence it makes sense to assign 
it some physical meaning. This is certainly consistent 
with previous HF studies using the ST0-3G minimal ba- 
sis set-^- and the more extensive 6-3 IG** basis set^^I], 
Since our 6-311-1— |-G**(2d,2p) basis set is even larger, it 
is not too surprising that we find an additional unoccu- 
pied orbital, namely the 7ai(3s) orbital shown in Fig. |5] 
Although apparently at least partially localized, this or- 
bital remains unbound at all bond angles in the Walsh 
diagram and care should be taken not to overinterpret 
its physical nature. Nevertheless this 7ai(3s) orbital in- 
tervenes in an important way in the interpretation of our 
calculated electronic absorption spectra and it is upon 
analysis of the spectra that we will be able to associate 
this orbital with the 3s Rydberg state. 

One should also be conscious in using the calculated 
TDB3LYP absorption spectrum to assign the experimen- 
tal gas phase UV spectrum that the TDDFT ionization 
continuum begins at minus the value of the HOMO or- 
bital energy^. The value of — ehomo obtained in the 
different calculations is 6.40 eV with the LDA, 7.68 eV 
with the B3LYP hybrid functional, and 12.27 eV with 
HF. As expected from Koopmans' theorem, the HF value 
of — enoMO is in reasonable agreement with the experi- 
mental ionization potential of 10.57 eV^. The presence 
of a fraction of HF exchange in the B3LYP hybrid func- 
tional helps to explain why its value of — ehomo lies be- 
tween that of the pure DFT LDA and that of HF. As far 
as our TDB3LYP calculations are concerned, the value 
of — ehomo means that assignment of experimental ex- 
citation energies higher than 7.7 eV should be avoided. 
Fortunately this still allows us to assign the first singlet 
excitation energies. Examination of the TDB3LYP co- 
efficients yields the following assignments for the three 
principal UV absorption peaks: 

6.69 eV : 1^ Bi[2bi{n) ^ 7ai{3s)] (5.2) 
7.14 eV : 2^ Bi[2bi{n) ^ SaiiSp,)] (5.3) 
7.36 eV : 2^1 [26i(n) ^ 36i(3p^)] . (5.4) 

Comparison with the experimental assignment of Basch 
et al. justifies the identification of these orbitals with the 
Rydberg orbitals 3s, 3pz, and 3px- 

It is worth pointing out that the expected 
^[26i(n), 462(0'*)] excitation is of ^^2 symmetry 
and as such corresponds to a spectroscopically forbidden 
transition. It is in any event fairly high in energy (9.77 
eV with TDHF). The ^[6ai(cr), 462(0"*)] transition has 
symmetry and so is spectroscopically allowed, but is 
still found at fairly high energy in our calculations (8.15 
eV with TDLDA, 8.37 with TDB3LYP, and 9.93 with 
TDHF). 



B. C2v Ring-Opening 

We now consider how DFT performs for describing the 
ground state and how well TDDFT performs for describ- 
ing the lowest excited state of each symmetry for 
ring-opening of oxirane. Comparisons are made against 
CASSCF and against high-quahty DMC energies calcu- 
lated at the same geometries. All geometries along the 
pathway have been fully optimized at each O - C - O 
ring-opening angle (a) using the B3LYP functional. The 
orbital energies as a function of ring opening angle have 
already been given in the Walsh diagram (Fig. |4]). 

Fig. [6] shows the TDB3LYP and TDLDA curves for 
the ground (l^^i) state and the lowest excited-states of 
each symmetry (21^1, l^Ai, l^Bi, l^Bi, l^^a, 1^^2, 
11^2, and 1^B2). Several things are worth noting here. 
The first point is that, some of the differences between 
the TDB3LYP and TDLDA curves are apparent, not 
real, as the points for the two graphs were not calculated 
at exactly the same angles. In particular, the TDLDA 
misses the point at 120° where we encountered serious 
convergence difficulties due to a quasidegeneracy of a 
and (7* orbitals (Fig.|4|. Under these circumstances, the 
HOMO, which suffers from self- interaction errors, can lie 
higher than the self- interaction- free LUMO. As the pro- 
gram tries to fill the orbitals according to the usual auf- 
bau principle, electron density sloshes on each iteration 
between the two orbitals making convergence impossible 
without special algorithms. The TDB3LYP calculations 
were found to be easier to converge, presumably because 
they have less self-interaction error. 

The second point is that only deMon2k TDLDA cal- 
culations are reported here, although we have also calcu- 
lated TDLDA curves using Gaussian and found similar 
results when both programs printed out the same infor- 
mation. Unfortunately, GAUSSIAN did not always print 
out the lowest triplet excitation energy, so we prefer to 
report our more complete deMon2k results. 

The third point is the presence of a triplet instability. 
This means that, going away from the equilibrium ge- 
ometry, the square of the first triplet excitation energy 
decreases, becomes zero, and then negative. The exact 
meaning of triplet instabilities will be discussed further 
below. For now, note that we follow the usual practice 
for response calculations by indicating an imaginary ex- 
citation energy as a negative excitation energy. However 
it is important to keep in mind that a negative excita- 
tion energy in this context is only a common convention 
and not a physical reality. Incidentally, we presume that 
triplet instability is the source of the difficulty with the 
Gaussian output mentioned above-the associated coef- 
ficients become more complicated and Gaussian may 
have difficulty analyzing them. 

A fourth and final point is that the TDDFT ionization 
threshold occurs at minus the value of the HOMO en- 
ergy which is significantly underestimated with ordinary 
functionals such as the LDA and B3LYP^^. This arti- 
ficially low ionization threshold is indicated in the two 
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FIG. 6: C2v ring opening curves: ground state (l^^i) curve calculated using the B3LYP (Gaussian) or the LDA (deMon2k) 
functional, lowest excited state curve of each symmetry {2^A\, I'^^i, l^Bi, I'^Bi, 1^^2, l'^^2, 1^^-82, and l^i?2) calculated 
using the TDB3LYP (Gaussian) and the TDLDA (deMon2k) excitation energies added to the B3LYP and the LDA ground 
state energy, respectively. The energy zero has been chosen to be the ground state energy for the 60° structure. Note that the 
"negative excitation energies" for the l'^i32 state relative to the ground state are really imaginary excitation energies (see text). 
Also shown is the TDDFT ionization threshold at — enoMO- 



graphs in Fig. |6] In the TDB3LYP case, the TDDFT 
ionization threshold is high enough that it is not a par- 
ticular worry. In the TDLDA case, the TDDFT ioniza- 
tion threshold is lower by about 1 eV. This may explain 
some of the quantitative differences between the high- 
lying TDB3LYP and TDLDA curves, although by and 
large the two calculations give results in reasonable qual- 
itative, and even semi-quantitative, agreement. 

We now wish to interpret the TDDFT curves. TDDFT 
excitations may be characterized in terms of single elec- 
tron excitations from occupied to unoccupied MOs^. 
Unlike in HF, unoccupied and occupied orbitals in pure 
DFT (e.g., the LDA) see very similar potentials and 
hence the same number of electrons. This means that 
the orbitals are preprepared for describing electron ex- 
citations and that simple orbital energy differences are 
often a good first approximation to describing excitation 
energies. In the two-orbital model and pure DFT, the 
singlet, ujs, and triplet, lot, TDDFT excitation energies 
for the transition from orbital i to orbital a are. 



integrals that. 



LOS 



LOT 



FT.il 



(ea-e») {ea-e,) + 2{ia\fj - f^-^\ai) .(5.5) 



In the TDA. this becomes 



(^s)tda = i^a - e{) + {ia\2fH + /i^' + /ic^l"*) 
(^t)tda ^ i^a - £») + {ia\flc - flc\ai) , (5.6) 

from which it is clear from the sizes and signs of the 



('^t)tda < - < i^s) 



TDA 



(5.7) 



For Rydberg states the inequalities become near equali- 
ties because the overlap of orbitals i and a goes to zero. 
While no longer strictly valid, these general ideas are still 
good starting points for understanding results obtained 
with the hybrid functional B3LYP. 

The frontier MOs shown in the Walsh diagram (Fig.|4]) 
supplemented with the addition of the two additional Ry- 
dberg orbitals, 3bi{3px) and 5^2 (3py), provides a use- 
ful model not only for interpreting the results of our 
TDDFT calculations, but also for constructing the ac- 
tive space necessary for the CASSCF calculations used 
to construct the QMC wavef unctions. Fig. [7] shows the 
results of our best QMC calculation and of the results 
of the CAS(4,6) calculation on which it is based. For 
the most part, the CAS and DMC curves differ quan- 
titatively but not qualitatively. The exceptions are the 
I'^Ai and l^i?2 curves where the DMC results, which are 
significantly lower than the corresponding CAS(4,6) re- 
sults at some geometries, contain important amounts of 
dynamic correlation not present at the CAS(4,6) level of 
calculation. In the remaining graphs (except for Fig. [9]) 
we will suppress CAS (4,6) curves in favor of presenting 
only the higher quality DMC curves since these offer the 
better comparison with TDDFT. 

We now give a detailed comparison of our TDDFT re- 
sults with DMC methods, dividing the states into three 
sets. The first set consists of states where double excita- 
tions (lacking in the TDDFT adiabatic approximation) 
are likely to be important. The second set consists of 
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FIG. 7: C2v ring opening curves calculated with CASSCF and DMC: curves for the lowest state of each symmetry (l^Ai, 
l^^i, l^Bi, l^Bi, 1M2, l^A2, 1^-62, and 1^52) are calculated using CAS(4,6) without state averaging, while the 2^ Ai is the 
result of adding the excitation energy from a state averaged calculation to the ground state l^Ai curve calculated without 
state averaging. The energy zero has been chosen to be the ground state energy for the 60° structure. Note that the negative 
excitation energies for the l^i32 state relative to the ground state are really negative excitation energies. Numerical DMC 
energies are listed in Appendix [B] 



states which show the effects of triplet instabilities. And 
the third set consists of Rydberg excitations. 

The states of Ai symmetry are those most likely to be 
affected by the absence of double excitations. Although 
not shown here, we have constructed the two HF curves 
obtained by enforcing a double occupancy of either the 
&ai{a) or the 462(0-*) orbital. Both single determinants 
have ^Ai symmetry. The two curves thus generated sim- 
ply cross at about 120°. In the absence of configuration 
mixing, the ground state curve always follows the lower 
state and shows an important cusp at 120°. Introducing 
configuration mixing via a CASSCF calculation leads to 
an avoided crossing. 

Now, DFT is different from HF because DFT is exact 
when the functional Exc is exact while HF always remains 
an approximation. On the other hand, the Kohn-Sham 
equations of modern DFT resemble the HF equations and 
tend to inherit some of their faults when approximate 
functionals are used. Fig.[8]shows the comparable curves 
at the TDDFT and DMC levels. As expected, the un- 
physical cusp is absent in the l^Ai ground state curve 
at the DMC level of calculation. The unphysical cusp is 
present at both the B3LYP and LDA levels where sig- 
nificant divergences between the DFT and DMC calcula- 
tions occur between about 100° and 150°. However, this 
region is very much reduced compared to what we have 
observed to happen for the HF curve where significant 
divergences beginning at about 75°. This is consistent 
with the idea that even DFT with approximate function- 
als still includes a large degree of electron correlation. 

As to the 2^Ai excited state curve, only the DMC 
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FIG. 8: C21, ring opening curves: 1^ Ai and 2^ Ai states. Note 
that the 1^ Ai state has also been included at the DMC level 
of calculation. 



calculation shows an indication of an avoided crossing. 
The inclusion of the I'^^i excited state curve calculated 
at the DMC level helps to give a more complete under- 
standing of the curve for this state. Below about 100° 
and above about 150°, the 2^Ai and l^Ai states have 
nearly the same energy, consistent with the idea that 
these correspond to the ^[ai(cr), ai(3s)] Rydberg tran- 
sition below 110° and to the ^ [62(0-*), 62(8^2)] Rydberg 
transition above 150°. In between the ^[ai(a)^, 62(0'*)^] 
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FIG. 9: C2v ring opening curves: 1^B2 and 1^B2 states. Note 
that the ground state (l^Ai) curve has only been shown for 
the B3LYP calculation since the LDA curve is nearly identical. 
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FIG. 10: C2v ring opening curves: l^Bi state. The TDLDA 
and TDLDA/TDA curves are practically superimposed. 



two-electron valence excitation cuts across the ^Ai man- 
ifold. The TDDFT calculations are qualitatively capable 
of describing these one-electron Rydberg excitations but 
not of describing the two-electron excitation in the bond- 
breaking region. As expected, the cusp in the DFT l^Ai 
ground state curve is simply reflected as cusps in the 
TDDFT excited-state curve. 

The only state which shows triplet instabilities is the 
1^B2 state. The TDLDA 1^B2 energies are imaginary 
(even if designated as negative in Figs. |6] and [9| between 
about about 100° and about 140° while the TDB3LYP 
1^B2 energies are imaginary over a larger range, between 
about 90° and 160°. Certainly, one way to under stand 
how this can happen is through the formulae Eq. (5.51, 



but a better way to understand triplet instabilities is in 
terms of the fact that the quality of response theory en- 
ergies depends upon having a high-quality ground state. 
Problems occur in TDDFT excitation energies because 
the functional E^c is only approximate. 

Stability analysis and triplet instabilities have already 
been discussed in a general way in Sec. |IT] where it was 
seen that no symmetry breaking is expected for a closed- 
shell singlet ground state when E^c is exact. However, 
symmetry breaking can occur when E^^ is approximate. 
To deepen our understanding of triplet instabilities, let 
us look at a two-orbital model for the dissociation of the 
a bond in H2, which is similar to the dissociation of the 
CC a bond in oxirane. Here we follow the argument in 
Ref.'22l and consider the wave function. 



iVl - A2cr + Act*, y^T^^a -~Xa*\. 



(5.8) 



which becomes |syi, sb| as A ^ l/\/2. The corresponding 
energy expression is. 



Ex — En 



2y 



0{}?). 



(5.9) 



This result suggests the more general result—® that sym- 
metry breaking will occur if and only if there is an imag- 
inary triplet excitation energy {u)\ < 0). Certainly, one 
way to overcome the problem of triplet excitation ener- 
gies is to improve the quality of the exchange-correlation 
functional. However, another way is to use the TDA. 
The reason is clear in HF theory where TDHF/TDA is 
the same as CIS, whose excited states are a rigorous up- 
per bound to the true excitation energies. Therefore, not 
only will the excitations remain real but variational col- 
lapse will not occur. To our knowledge, there is no way 
to extend these ideas to justify the use of the TDA in 
the context of TDDFT calculations except by carrying 
out explicit calculations to show that the TDA yields 
improved PESs for TDDFT calculations. This was pre- 
viously shown for H2'^^ and is evidently also true for oxi- 
rane judging from Fig. [9] where the TDLDA/TDA curve 
is remarkably similar to the DMC curve. 

As shown in Fig. [9] triplet instabilities are often also 
associated with singlet near instabilities. In this case, 
the TDDFT 1^B2 singlet excitation energies are much 
too low. The TDA brings the TDLDA curve into the 
same energetic region of space as the corresponding 
DMC curve, but still compared with the DMC curve the 
TDLDA TDA curve appears to be qualitatively incorrect 
after 120° where it is seen to be decreasing, rather than 
increasing, in energy as the angle opens. This, in fact, 
is the behavior observed in the CAS(4,6) singlet calcula- 
tion but is opposite to what happens in the more accurate 
DMC calculation. 

The remaining 1^'^Bi and 1^'^A2 states are primar- 
ily Rydberg in nature with the corresponding singlets 
and triplets being energetically nearly degenerate. The 
graph for a single one of these states suffices to illus- 
trate the general trend observed in the case of all four 
states. Fig. 10 shows what happens for the l^Bi states. 
All the curves have qualitatively the same form. This is 
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FIG. 11: Comparison of TDLDA/TDA and DMC C2v ring opening curves. 



particularly true for the TDB3LYP curve which closely 
resembles the DMC curve but shifted down by a little 
over one eV in energy. Notice also that the TDLDA and 
TDLDA/TDA curves are essentially superimposed. This 
is because, 



Thermal Reaction 
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(5.10) 



For Rydberg states the overlap between orbitals i and a 
is very small and so the difference between full and TDA 
TDLDA calculations can be neglected. From Eq. (5.4) 



this also means that the Rydberg excitation energies also 
reduce to simple orbital energy differences. 

shows that the TDLDA/TDA yields, at least 
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in this mostly qualitatively reasonable description 

of photochemically important PESs as compared with 
high-quality QMC PESs. By this we only mean that 
corresponding TDLDA/TDA and QMC potential energy 
curves have typically the same overall form and are found 
at roughly similar energies. As discussed above, signifi- 
cant differences do remain in the shapes of some potential 
energy curves {2^Ai and 1^-82) and the detailed energy 
ordering of the curves differs between the TDLDA/TDA 
and QMC results. 



Fig. 13 shows the ground state (^o) PES. In accor- 



dance with the WH model there appears to be a much 
lower energy barrier for con ring opening than for dis 
ring opening. Also shown is the lowest triplet state (Ti). 
It is clear that triplet instabilities occur along the dis- 
rotatory pathway as the CC a bond breaks. Remark- 
ably, they do not occur along the conrotatory pathway 
in the LDA. (They should, of course, be absent when 
the exchange-correlation functional is exact.) Although 
triplet instabilities account for about 50% of the LDA 
surface, this fraction actually increases to 93% for the 
B3LYP functional where triplet instabilities occur along 
both the conrotatory and disrotatory reaction pathways. 
Finally, when the HF method is used triplet instabilities 
account for 100% of the PES. Perhaps because we forced 
the OCII2 to lie in a single plane, triplet instabilities 
appear to be nearly everywhere in these last methods, 
potentially spelling trouble for mixed TDDFT/classical 
photodynamics calculations. We would thus like to cau- 
tion against the use of hybrid functionals for this type 
of application. At the same time, we reiterate our rec- 
ommendation to use the TDA because, of course, the 
TDLDA/TDA excited state PESs (Fig. K^ show no 



C. Con- and Disrotatory Ring Opening 

Conrotatory and disrotatory potential energy surfaces 
were calculated using the simplifying assumption that 
each set of 4 atoms OCH2 is constrained to lie in a plane. 
Our coordinate system is defined in Fig. 12 All other 



geometrical parameters were relaxed for the thermal re- 
action preserving Cg and C2 symmetry for the con and 
the dis path, respectively. 




FIG. 12: (a, 9) coordinate system used in this paper. 
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FIG. 13: TDLDA So (light grey) and Ti (dark grey) PESs. 
Note that a "negative excitation energy" is just a convenient 
graphical trick for representing an imaginary excitation en- 
ergy. Negative excitation energies correspond to triplet insta- 
bilities. 
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FIG. 14: TDLDA/TDA So (light grey) and Ti (dark grey) 
PESs. In this case negative excitation energies are real, not 
imaginary, quantities. 

triplet instabilities. 

2. Photochemical Reaction 

In considering whether the WH model has any validity 
for describing con- versus disrotatory rotation for pho- 
tochemical ring opening in oxirane, we are immediately 
faced with the problem that there are a large number of 
excited states with similar energies which cross each other 
(Fig. |7|. Provisionally, we have decided to assume that 
state symmetry is conserved and so to look only at one 
PES, namely the l^Bi^i surface (which we shall simply 
refer to as Si) which begins as the ^[26i(n) 7ai(3s)] 
excitation for the closed cycle and then evolves as the 
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FIG. 15: TDLDA/TDA So (light grey) and Si (white) poten- 
tial energy surfaces. 



molecular geometry changes. In so doing, we are follow- 
ing the spirit of simple WH theory which makes max- 
imum use of symmetry. Turning back to this theory 
(Fig. [T]) and thinking of excitation to a Rydberg state 
as analogous to electron removal, it suffices to remove an 
electron from the n orbital of the thermal WH diagram, 
thereby giving the u^n} configuration, to have an idea of 
what might be the relative importance of the con and dis 
mechanisms for the surface. In this case the con mech- 
anism leads to the energetically nearest-neighbor single 
excitation, ct — > n, and the dis mechanism leads to the 
nearest-neighbor single excitation, a ^ a* . Since both 
mechanisms correspond to nearest-neighbor excitation, 
no particular preference is a priori expected for one mech- 
anism over the other. The ground (6*0) and excited state 
singlet {Si) surfaces calculated with the TDLDA/TDA 
are shown in Fig. |15[ Indeed the simple WH theory ap- 
pears to be confirmed in that the Si surface is remarkably 
flat. 

Unfortunately this analysis is far too simplistic. In par- 
ticular. Kasha's rul^SHteHs us that the first excited triplet 
(Ti) or singlet {Si) states are the most hkely candidates 
for the initiation of a photochemical reaction, where now 
Ti and 5'i are the globally lowest states and not necessar- 
ily the lowest state of a given symmetry. This is based 
upon the idea that relaxation of higher excited states is 
rapid. Such relaxation is due to environmental effects or 
vibronic coupling which need not preserve the symme- 
try of the electronic state. Hence at one geometry is 
not necessarily Si at another geometry. Indeed looking 
again at the DMC curves in Fig. [7] it is easy to believe 
that the molecule will eventually arrive in the 1^B2{<J, a*) 
state during ring opening and that, because the a orbital 
is higher in energy than the a* orbital at this geometry 
(Fig. [4]), that the usual WH argument will still predict a 
preference for the dis mechanism. 

Again we are falling into a trap imposed by the use of 
symmetry. Excitation, for example, from the 26i(n) or- 
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bital (Fig.[3| into the 8ai(3p^) Rydberg orbital (Fig. [s]) 
might lead to preferential CO, rather than CC, bond 
breaking to the extent that ring opening augments the 
valence CO (cr*) character of the target orbitaP^. The 
difficulty of predicting a priori such behavior is part of 
what motivates us to move towards dynamics as a better 
tool for photochemical modeling. 



VI. CONCLUSION 

In this paper, we have examined the potential energy 
curves and surfaces for the symmetric ring opening of 
oxirane to assess possible difficulties which might be en- 
countered when using TDDFT in photodynamics simu- 
lations of its photochemistry. Our TDDFT calculations 
provided useful insight helpful in constructing the ac- 
tive space needed for more accurate CASSCF and still 
more accurate DMC calculations. It is indeed worth not- 
ing that identifying active spaces is probably one of the 
more common uses of TDDFT in photochemical studies. 
Here, we summarize our main conclusions obtained by 
comparing our TDDFT and DMC results. 

Oxirane does not seem to be a molecule where charge 
transfer excitations are important. The artificially low 
ionization threshold typical in TDDFT for most func- 
tional is still high enough for the B3LYP functional so 
as not to pose a serious problem. Even for the LDA func- 
tional, where the TDDFT ionization threshold is lower, 
the shapes of the excited-state Rydberg curves seem to 
be qualitatively correct even if they should not be con- 
sidered quantitative. 

Problems do show up as the CC a bond breaks. By 
far, the most severe problem encountered is the presence 
of triplet instabilities and singlet near instabilities, where 
the excited-state PES takes an unphysical dive in energy 
towards the ground state. In the case of the triplet, the 
excitation energy may even become imaginary, indicating 
that the ground state energy could be further lowered by 
allowing the Kohn-Sham orbitals to break symmetry. Al- 
though this problem is very much diminished compared 
to that seen in the HF ground state, it is still very much 
a problem as seen for example in the 2D surfaces of the 
con and disrotatory ring opening of oxirane. 

It is difficult to overstate the gravity of the triplet in- 
stability problem for the use of TDDFT in photodynam- 
ics simulations. One should not allow symmetry break- 
ing as it is an artifact which should not occur when the 
exchange-correlation functional is exact. Moreover, there 
may be more than one way to lower the molecular energy 
by breaking symmetry, and searching for the lowest en- 
ergy symmetry broken solution can be time consuming. 
Finally, the assignment of excited states using broken 
symmetry orbitals is far from evident. On the whole, it 
thus seems better to avoid the problem by using the TDA 
to decouple (at least partially) the quality of the excited- 
state PES from that of the ground state PES. Our cal- 
culations show that this works remarkably well for the 



l^i?2 curve along the C2v pathway, in comparison with 
good DMC results. The l^i?2 state, which appeared to 
collapse in energy without the TDA in the region of bond 
breaking, is also restored to a more reasonable range of 
energies. For this reason, we cannot even imagine carry- 
ing out photochemical simulations with TDDFT without 
the use of the TDA. Or, to put it more bluntly, the TDA, 
which is often regarded as an approximation on conven- 
tional TDDFT calculations, gives better results than does 
conventional TDDFT when it comes to excited-states 
PESs in situations where bond breaking occurs. 

The TDA is unable to solve another problem which oc- 
curs as the CC a bond breaks, namely the presence of an 
unphysical cusp on the ground state curve (or surface). 
This cusp is there because of the difficulty of approximate 
density functionals to describe a biradical structure with 
a single determinantal wave function. The cusp is also 
translated up onto the excited-state curves because of the 
way that PESs are constructed in TDDFT. However the 
problem is very much diminished compared to that seen 
in HF because of the presence of some correlation in DFT 
even when the exchange-correlation functional is approx- 
imate. The ultimate solution to the cusp problem is most 
likely some sort of explicit incorporation into TDDFT of 
two- and higher-electron excitations. Among the meth- 
ods proposed for doing just this are dressed TDDFT or 
polariza tion pr opagator corrections"^^ ^"^ and spin-flip 

Our examination of the con- and disrotatory ring- 
opening pathways revealed another important point 
which has nothing to do with TDDFT. This is that the 
manifold of excited state PESs is too complicated to in- 
terpret easily. The best way around this is to move from 
a PES interpretation to a pathway interpretation of pho- 
tochemistry, and the best way to find the pathways mov- 
ing along and between adiabatic PESs seems to be to 
carry out dynamics. For now, it looks as though the 
use of TDDFT/TDA in photodynamics simulations may 
suffice to describe the principal photochemical processes 
in oxirane. On-going calculation J^Sl do indeed seem to 
confirm this assertion in so far as trajectories have been 
found in g ood ag reement with the accepted Gomer-Noyes 
mcchanisrrPlESl, 



APPENDIX A: PHOTOCHEMISTRY OF 
OXIRANES 

The generic structure of oxiranes is shown as the start- 
ing point in the chemical reactions in Figs. [16] and |17[ 
When Ri, R2, R3, and R4 are hydrogens or alkyl groups, 
then the prefered reaction is CO cleavage both photo- 
chemically and thermally (Step 1, Fig. 16). In particular 
it is estimated that the CO rupture energy is about 52 
kcal/mol while the CC rupture energy is 5-7 kcal/mol 
higher Since the molecule is not symmetric along the 
CO ring-opening pathway, the WH model does not ap- 
ply. Photochemical CO ring opening may be followed 
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TABLE II: Ai DMC energies as a function of COC ring- 
opening angle along the C2v ring-opening pathway. 



Ai DMC Energies and statistical error in Hartree 



Z COC (°) 


I Ai 


VAi 


60. 


0.0000 (0.0009) 


0.3099 (0.0010) 


75. 


0.0212 (0.0009) 


0.2811 (0.0010) 


90. 


0.0651 (0.0010) 


0.2836 (0.0009) 


105. 


0.0998 (0.0010) 


0.2998 (0.0010) 


120. 


0.1222 (0.0009) 


0.3406 (0.0011) 


135. 


0.1153 (0.0009) 


0.3223 (0.0009) 


150. 


0.1083 (0.0009) 


0.3106 (0.0010) 


165. 


0.1048 (0.0010) 


0.3025 (0.0010) 


179.5 


0.1017 (0.0010) 


0.3007 (0.0010) 



FIG. 16: Typical reactions of alkyl oxiranes. 



o (+) 

R, \ 

R4 



turns out to be not straightforward, particularly in the 
photochemical case. There are several reasons for this. 
First of all, too much asymmetry should be avoided in or- 
der to assure the applicability of the WH orbital symme- 
try conservation rule and so avoid passing directly onto 
carbene formation (Step 2, Fig. 17 1. Secondly, the sub- 



stituted oxirane should include groups which are bulky 
enough to confer the structural rigidity needed to avoid 
premature radiationless relaxation, but not bulky enough 
to favor carbene formation. The predictions of the WH 
model have been found to hold for cis- and trans-1,2- 
diphenyloxiranc'''-^ w hile carbene formation dominates 
for tetraphenyl oxiran^^^. 



APPENDIX B: SUPPLEMENTARY MATERIAL 



FIG. 17: Typical reactions of aryl oxiranes. 



by alkyl migratioEp^ (Step 2, Fig. 16 1. In the particu- 
lar case of oxirane itself (Ri=R2=R3=R4=H), hydrogen 
migration is foll owed by breaking of the CC single bond 
(Step 3, Fig. 16). This is the Gomer-Noyes mechanismSSj 
which was confirmed experimentally by Ibuki, Inasaki, 
and TakesakP. 

In contrast, cyano and aryl substitutions favor CC 
bond breaking (Step 1, Fig. 17) to form what is often 
refered to as a 1,3-dipolar species or a carbonyl yhde. 
This is the case where there may be sufficient symme- 
try that the WH model applies. The photochemistry 
of phenyl and phenyl substituted oxiranes has been re- 
viewed in Refs) ^6 | 77 | 78 | j^jjj pj^ggg 565-566 of Rei]^. 

Evidence for carbonyl ylides goes back to at least the 
1960s when UUman and Milks invest igated t he tautomer- 
ization of 2,3-diphenylindenone oxide^^ZHMl ^nd Linn and 
Benson investigate d the ri ng-opening reaction of tetra- 
cyanoethylene oxid^^^llinil Finding cases where the WH 
model applies and where its predictions can be verified 



Benchmark quality diffusion Monte Carlo (DMC) en 
ergies calculated at the geometries given in Table IVH 
are reported here for the ground state {l^Ai, Tables |II 
and HI I and the lowest excited state of each symm etry 
l^Ai, Table [h] l^Bi, l^Bi, Table |lV 



Table 



HI 



(2'^. _ - _ - - _ 

1^A2, 2^A2, Table|v] and l^Bs, l^-Ba, Table [Vl|. We set 

the zero of the energy to coincide with the DMC l^Ai 
energy at 60° and report the statistical error on the en- 
ergy in parenthesis. Note that the negative value of the 
l^Ai SA-DMC energy at 60° is statistically compatible 
with zero. We hope that these DMC data will encourage 
further developing and testing of improved TDDFT al- 
gorithms suitable for addressing the problems mentioned 
in this article. 
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TABLE III: Ml SA-DMC energies as a function of COC ring 
opening angle along the C2v ring-opening pathway. 

^ Ai SA-DMC Energies and statistical error in Hartree 



Z COC (°) 


1^1 


2^1 


60. 


-0.0018 (0.0009) 


0.3246 (0.0009) 


75. 


0.0254 (0.0010) 


0.2904 (0.0009) 


90. 


0.0694 (0.0010) 


0.2865 (0.0010) 


105. 


0.1011 (0.0010) 


0.2962 (0.0009) 


120. 


0.1263 (0.0009) 


0.2699 (0.0010) 


135. 


0.1196 (0.0010) 


0.2902 (0.0010) 


150. 


0.1103 (0.0010) 


0.3109 (0.0011) 


165. 


0.1054 (0.0010) 


0.3063 (0.0010) 


179.5 


0.1037 (0.0010) 


0.3054 (0.0010) 



TABLE IV: Bi DMC energies as a function of COC ring 
opening angle along the C2v ring-opening pathway. 

Bi DMC Energies and statistical error in Hartree 



Z COC (°) 


l^Bi 


l^Bi 


60. 


0.2788 (0.0010) 


0.2636 (0.0010) 


75. 


0.3083 (0.0009) 


0.3026 (0.0010) 


90. 


0.3138 (0.0010) 


0.3140 (0.0009) 


105. 


0.3329 (0.0010) 


0.3316 (0.0010) 


120. 


0.3058 (0.0010) 


0.2970 (0.0010) 


135. 


0.3023 (0.0010) 


0.2887 (0.0010) 


150. 


0.3113 (0.0010) 


0.2972 (0.0010) 


165. 


0.3274 (0.0010) 


0.3111 (0.0011) 


179.5 


0.3343 (0.0010) 


0.3155 (0.0010) 



from the Mexican Ministry of Education via a CONA- 
CYT (SFERE 2004) scholarship and from the Universi- 
dad de las Americas Puebla (UDLAP). 

We Hke to acknowledge useful discussions with Math- 
ieu Maurin, Enrico Tapavicza, Neepa Maitra, Ivano Tav- 
erneUi, Todd Martinez, and Massimo Ohvucci. 

MEC and FC thank Pierre Vatton, Denis Chara- 
poff, Regis Gras, Sebastien Morin, and Marie-Louise 



Dheu-Andries for technical support of the Departement 
de Chimie Molecularie and Centre d' Experimentation 
le Calcul Intensif en Chimie (CECIC) computers. CF 
aknowledges the support by the Stichting Nationale 
Computerfaciliteiten (NCF-NWO) for the use of the 
SARA supercomputer facilities. 



TABLE V: A2 DMC energies as a function of COC ring 
opening angle along the C2v ring-opening pathway. 

A2 DMC Energies and statistical error in Hartree 



Z COC (°) 


1^2 


1^2 


60. 


0.2972 (0.0010) 


0.3001 (0.0009) 


75. 


0.3021 (0.0010) 


0.2971 (0.0009) 


90. 


0.2901 (0.0010) 


0.2857 (0.0010) 


105. 


0.2835 (0.0010) 


0.2769 (0.0010) 


120. 


0.2870 (0.0010) 


0.2843 (0.0010) 


135. 


0.3344 (0.0010) 


0.3319 (0.0011) 


150. 


0.3179 (0.0010) 


0.3165 (0.0010) 


165. 


0.3111 (0.0010) 


0.3088 (0.0010) 


179.5 


0.3086 (0.0010) 


0.3066 (0.0010) 



TABLE VI: B2 DMC energies as a function of COC ring 
opening angle along the C2v ring-opening pathway. 

B2 DMC Energies and statistical error in Hartree 



Z COC 1 


n i'B2 


1^52 


60. 


0.3484 (0.0009) 


0.3301 (0.0010) 


75. 


0.3040 (0.0010) 


0.2395 (0.0010) 


90. 


0.2810 (0.0010) 


0.1657 (0.0010) 


105. 


0.2549 (0.0011) 


0.1229 (0.0009) 


120. 


0.2470 (0.0011) 


0.1165 (0.0010) 


135. 


0.2481 (0.0010) 


0.1232 (0.0009) 


150. 


0.2604 (0.0010) 


0.1410 (0.0010) 


165. 


0.2718 (0.0010) 


0.1611 (0.0009) 


179.5 


0.2733 (0.0010) 


0.1667 (0.0010) 
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